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Numerically efficient transfer matrix technique for studying statistics of coherent adsorbates on 
small nanotubes has been developed. In the framework of a realistic microscopic model fitted to 
the data of ab initio calculations taken from literature sources, the ordering of potassium adsorbate 
on (6,0) single-walled carbon nanotube has been studied. Special attention has been payed to the 
phase transition-like abrupt changes seen in the adsorption isotherms at low temperature. It has 
been found that the behavior during the transitions conforms with the universality hypothesis of the 
theory of critical phenomena and is qualitatively the same as in the one dimensional Ising model. 
Quantitatively the critical behavior can be fully described by two parameters. Their qualitative 
connection with the properties of interphase boundaries is suggested but further research is needed 
to develop a quantitative theory. 



PACS numbers: 64.70. Nd,68.43.-h 



I. INTRODUCTION 



A considerable part of the ongoing research on adsorp- 
tion in carbon nanostructures is driven by the problem 
of hydrogen storage at ambient conditions. 1-4 In particu- 
lar, the metallic adsorbates are expected to considerably 
enhance the hydrogen uptake^ - — because the storage ca- 
pacity of purely carbon structures is insufficient from a 
practical point of view 4 The adsorption of gases£~— al- 
lows, inter alia, to gain deeper insight into the depen- 
dence of sorption on various characteristics of adsorbate 
molecules, such as their sizeJ^- 

From the storage perspective, the most promising 
among carbon nanostructures are the single-walled nan- 
otubes (SWNTs) because of their large surface to weight 
ratio^ Since the storage capacity is defined mainly by 
the adsorbing surface^ theoretical studies of the ad- 
sorption for simplicity are often performed on individual 
SWNTs^-'ii'i 2 - The hydrogen uptake predicted in such 
studies is sometimes very high£i£ii2 but their significance 
for the storage is not clear because the calculations are 
usually made for periodic structures at zero temperature 
with only crude estimates of temperature effects some- 
times being m&de&H 

Temperature effects, however, may strongly influence 
predictions based on zero-temperature calculations. For 
example, at finite temperatures the ordered structures 
cannot exist in one-dimensional systems in the thermo- 
dynamic limiti 1 ^ Instead, if temperature is sufficiently 
low, a disordered state is formed with extended local or- 
der corresponding to the T = K ordered structure. 
From continuity considerations it is reasonable to assume 
that at sufficiently low temperature this quasi-ordered 
structure should be as good a hydrogen absorber as the 
zero-temperature one. Thus, from the storage point of 
view the question is how large are the temperatures at 
which the zero-temperature predictions can still be re- 
lied upon. In the closely related problem of adsorption 



on the two-dimensional (2D) surface this question can 
be answered with the help of phase diagrams where or- 
dered and disordered phases are separated by well defined 
boundaries j 14 ' 15 

The aim of the present paper is to adopt the techniques 
of Refs. [i3l5l where adsorption of hydrogen and oxygen 
on 2D surfaces were investigated to the case of adsorption 
on individual SWNTs and to find out what can be said 
about the quasi-ordered structures in the absence of well 
defined finitc-tcmperature phase boundaries. To imple- 
ment this approach, some assumptions and approxima- 
tions need be made which are usually specific to the type 
of the adsorbate under consideration. For concreteness, 
we will discuss them using as an example the potassium 
deposit on (6,0) zigzag SWNT. This choice was moti- 
vated mainly by the fact that this system was studied 
with an ab initio technique in Ref. [6j where the ground 
state energies of six periodic structures were calculated. 
Such information is necessary for the implementation of 
the cluster expansion method (CEM) 1 ^— which allows 
one to derive an effective lattice gas Hamiltonian to use 
in the solution of statistical problems. In connection with 
Refs. [l6T - fl8l which deal with binary alloys it is pertinent 
to point out that the lattice gas model is formally equiv- 
alent to the binary alloy which allows for the use of tech- 
niques developed in the alloy theory to coherent surface 
adsorbates. It should also be noted that in non-metallic 
systems the CEM can be developed on the basis of the 
energies calculated with the use of model potentials. 11 
Furthermore, the effective Hamiltonian can be derived 
via fit to experimental data j 14 ' 15 

An important assumption made in the application of 
CEM to surface structures is that the adsorbate is in 
registry with the substrate lattice. In reality, however, 
this depends on the relative strength of interactions of 
adsorbate atoms between themselves and with the sub- 
strate. If the latter interaction is weak (which is the 
case in the system under consideration 1 ^) the coherence 
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with the substrate in sufficiently dense structures can be 
lost3 

Another difficulty is due to the substantial coverage- 
dependent charge transfer which strongly influences the 
interactions of adsorbate atoms with the substrate and 
with each other (see detailed discussion for the potas- 
sium depos its on graphite, graphene, and on SWNTs in 
Refs. 1 1 9U 2 lh - In principle, charge transfer should be ade- 
quately accounted for in the ab initio calculations. 6 But 
the system under consideration is apparently rather sin- 
gular because the adsorption energy of potassium atom 
on the graphen varies in different calculations from 0.44 
to 2.0 eVi^i On that scale our neglect of phonons whose 
characteristic Debye energy is usually an order of magni- 
tude smaller looks justified. The question remains on the 
importance of the entropic contribution at finite temper- 
ature due to the atomic vibrations. In the case of alloys 
this problem was reviewed in Ref. [22| where it was con- 
cluded that in the majority of cases the classical approx- 
imation should be sufficient. At least, this is justifiable 
in the case of atoms with large atomic mass like potas- 
sium. The classical statistical averaging in the harmonic 
approximation on which the phonon theory is based re- 
duces to the Gaussian integration over atomic coordi- 
nates which can be performed exactly. As explained in 
Appendix B of Ref. [23|, the energy part of the configura- 
tion free energy thus obtained coincides with the energy 
minimum at zero temperature. Thus, the equilibrium 
atomic positions obtained in the ab initio calculations 6 
at zero temperature provide the necessary average energy 
while the entropic part can be unified with the interac- 
tion part as an effective temperature-dependent contribu- 
tion into the pair interatomic interaction ! i 23 Because in 
the present study we are going to consider temperatures 
which are small in comparison with the pair interactions, 
we will neglect this contribution in our calculations. 

Finally, when deposited atoms or molecules are very 
light (He and H2 being the most important examples), 
quantum corrections became important at low tempera- 
tures and quantum treatment is preferable . 24 ' 25 In Ref. 
[241 . however, it was shown for the hydrogen molecules 
adsorbed in the nanotube bundles that a fully classi- 
cal regime sets in already at temperatures above 60 K. 
Because for the storage purposes the temperatures be- 
low the liquid nitrogen boiling point (77 K) are of little 
interest; 4 - the quantum corrections can be neglected in 
the studies oriented on storage applications. 

Thus, the main difficulties in the statistical description 
of the adlayers on SWNTs are due to the incommensu- 
rate structures and the poor accuracy of the interaction 
parameters. The accuracy can be improved either with 
the use of a better ab initio approach or by a direct fit 
to experimental data.— The loss of coherence with the 
substrate is a more serious problem because the lattice 
gas formalism usually requires a regular lattice to exist 
(see, however, Sec. 5.3 of Ref. [iH) . In adsorbates this usu- 
ally restricts the coverages at which the system retains 
its coherence to low values < 0.54^ 



But on the other hand, as is well known (see, e. g., 
Ref. Il7l ) , the lattice gas model is equivalent to the Ising 
model which is famous for being capable of describ- 
ing such disparate critical phenomena as the magnetic 
ordering in uniaxial magnets and the liquid-gas phase 
transition^ - — This similarity between the critical phe- 
nomena has been conceptualized in the universality hy- 
pothesis which has been amply confirmed by both exper- 
imental data and theoretical calculations in 2- and 3D 
systemsi 29 i 30 Hopefully it will work equally well in ID 
system s 31 ' 32 which may allow for the extension of our 
results to the incommensurate cases as well. 

In view of the many approximations and assumptions 
which need be accepted in order to implement the sta- 
tistical approach to the adsorption on SWNTs, in the 
present paper we will focus on the low temperature re- 
gions in the vicinity of critical points where the universal 
behavior sets in and where even significant inaccuracies 
in the microscopic description in most cases may be ir- 
relevant. 

In the next section we will explain the universality hy- 
pothesis for one-dimensional systems belonging to the 
Ising universality class; in Sec. Mil the effective Hamil- 
tonian will be derived in the framework of the CEM; in 
SeclIVIthe partition function will be calculated with the 
use of numerically efficient transfer matrix (TM) tech- 
nique and in Sec. [V] we present our conclusions. 

II. UNIVERSALITY IN ID 

Universality hypothesis constitutes one of pillars of the 
modern theory of critical phenomena^ - — It states that 
the singular part of the equations of state of all sys- 
tems belonging to the same universality class has the 
same functional form in the vicinity of the critical point. 
Unique for each system are only two constant parameters 
which define the scales of variation of the (dimensionless) 
external field 

L = h/k B T (1) 

and of the reduced temperature 

t = (T - T c )/T c , (2) 

where T c is the critical temperature. Thus, one may sim- 
plify the task of predicting the behavior of a system in the 
critical region by solving the simplest model belonging to 
the universality class of interest. The two parameters can 
be either found in independent calculations or derived 
from experimental data. For the general discussion of 
the universality we refer the interested reader to the vast 
literature on the subject^ - — while below we will consider 
only the Ising universality class in lDi^r— The peculiar- 
ity of this case is that there is no finite-temperature phase 
transitions in ID. Therefore, the approach to universality 
based on scaling variables (TTJ) and @ cannot be applied 
straightforwardly because T c = and the scaling variable 
t is undefined. 
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A solution to this problem was found in Ref. m. it 
was noted that instead of the scaling parameter ^ the 
correlation length £ can be used due to the relation 

t~C 1/l/ : (3) 

where v is the critical index which defines the divergence 
of the correlation length as £ ~ t~ v . Taking into account 
that in ID all critical indices are known exactly, on the 
basis of Eqs. (3.40)-(3.41) of Ref. l32l the equation of state 
in the scaling region can be written as 

M « Wm/2), (4) 

where M is the magnetization normalized as 

M(±oo) = ±1 (5) 

and W is the scaling function. The latter is universal 
for all systems belonging to the Ising universality class 
in ID except for two constant factors: one factor multi- 
plying W thus changing the range of variation of M in 
Eq. (|5J| and another one before its argument. We used 
this arbitrariness in Eq. (j4| by dividing the argument by 
two in comparison with Ref. 32. This definition is more 
appropriate to our purposes and according to Eq. ([3]) it 
does not change the scaling relations which are invariant 
under rescalings. 

The above formalism can be easily refashioned to de- 
scribe the lattice gas model via the equivalence transfor- 
mation 

o-j = 2rii - 1, (6) 

where cj{ = ±1 is the Ising spin on site i and rii = 0, 1 the 
corresponding occupation number. From this identity it 
follows that 

hai = firii - p/2, (7) 

where p = 2h. 

The coverage is defined as the lattice gas density 

P = {m), (8) 

where the angular brackets denote statistical averaging 
and the dependence of p on the lattice site is absent be- 
cause the system is assumed to be homogeneous and the 
spontaneous symmetry breaking is absent in ID. 13 
According to Eqs. © and ([7J Eq. (QJ takes the form 

(p - p c )/(Ap/2) « W[(p - p c )i/k B T], (9) 

where p c is the chemical potential at the critical point 
and Ap = p + — p_ is the total change of the density 
during the transition. Because at the critical point the 
critical density p c ~ (p + + P-)/2, the left hand side of 
Eq. (0) varies in the same range as in Eq. ([5]). We note 
that to achieve this we had to divide p — p c on the right 
hand side of Eq. ^ by Ap/2, not by p c as suggested 



in Ref. [29]. In this way we fix one of the arbitrary scale 
factors in our problem. 

Thus, according to the universality principle the be- 
havior of the system near any critical point can be de- 
scribe with the use of only two constant scale factors 
provided the universal function W is known. The lat- 
ter can be calculated for the simplest possible model, the 
Ising model with the fist neighbor interactions being the 
most obvious choice. 



A. ID lattice gas model 

Exact solutions of the ID Ising model can be found 
in many places, for example, in Eq. (3.39) of Ref. [32| 
or in Ref. |34|. But below for completeness we present 
the solution of the equivalent lattice gas model with the 
use of a variant of the nonsymmetric transfer matrix 
techniqu o 35 ' 36 which in Sec. HVI will be generalized to the 
case of nanotubes. 

In the process of adsorption the number of atoms on 
the surface is governed by the chemical potential p which 
may be controlled by the gas pressure if the adsorption 
from gaseous phase takes place (see, e. g., Ref. U3) or by 
the concentration of the adsorbate in the solution in the 
case of adsorption from a liquid. Therefore, the natural 
choice is the grand ensemble formalism with the partition 
function 

5 = Tre-f" 1 , (10) 

rii— 0,1 

where j3 = 1/ksT and H the Hamiltonian; for brevity the 
term with the chemical potential —p'Yln n-i is considered 
to be included into H. u With the use of Eq. dTDJ) the 
coverage can be found as [cf. Eq. ©] 

p={PN)- 1 d\n~,/dp, (11) 

where N is the number of deposition sites. In the case 
of ID lattice gas with only nearest neighbor interaction 
V\ (which we assume to be attractive) Eq. (flO|) can be 
written as 

N-l 

E 1D = Tr e^ nN T] exp(f3p ni - (3V ini n l+1 ). (12) 

71,-0,1 A A 

»=1 

We assume free boundary conditions corresponding to 
nanotubes with open ends. According to Eq. (fTIZ|) . Em 
can be calculated via the N — 1-st power of the transfer 
matrix 

^ = ( 6 /3m ePbt-Vt) J ■ ( 13 ) 

In the thermodynamic limit the reduced free energy Eq. 
(1291) of the ID lattice gas is 

<t> 1D = -p- 1 \n\+, (14) 
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where A+ is the largest eigenvalue of f. The logarithm 
in Eq. (fT4"| can be cast into the form 



In A + = S + ln(cosh S + Vsinh 2 6 + e^), 



(15) 



where 



S = p(p-V 1 )/2 = P(p-p c )/2. (16) 
The coverage can be found as 



1 1 



p = j3 1 d\n\ + /dp=- + 



sinh S 



2 2 y/ sm h 2 (5 



(17) 



As is easy to see, at low temperature e@ Vl — > and p(S) 
tends to the ^-function. In other words, all variation of 
p{5) is restricted to a narrow interval of p ~ p c . That is 
why this region is so important. The interaction poten- 
tial and the temperature may vary in very broad ranges 
but the values p = to the left of the interval and p = 1 
to the right of it will remain the same. In other words, 
these ranges of variation of p do not provide much use- 
ful information on the microscopies of the system. In 
contrast, in the vicinity of p c the slope 



dp _ e^coahS 
dpi ~ 4 (sinh 2 6 + e^f/ 2 



(18) 



will vary very strongly with temperature, with V±, with p, 
etc., so all quantities of interest are most easily measured 
near this quasi-transition point. 

At p — p c p in Eq. (|17p is equal to 0.5. This means 
that both phases — p = and p = 1 — are present in the 
system in equal proportion. At low temperature accord- 
ing to Eqs. (j3"Tj) and (|T5|) the atoms are correlated at long 
distances so the system looks as an intermittent mixture 
of the pure phases separated by interphase boundaries 
(IPBs). In the model under consideration the boundary 
energy at T — is easy to calculate. In the Ising spin 
representation Eqs. (JB])-© the spin-spin interaction at 
h = (p = 0.5) is 

0V4)X^+i = VtY^nirii+x -piN/2. (19) 



At zero temperature the IPB will separate the region of 
spins up from the spin down region. According to Eq. 
(fT!)]) . in comparison with the ordered system the energy 
cost is 



E h 



\Vi\/2. 



(20) 



Furthermore, in a system with free boundaries there are 
two equally probable possibilities for an IPB: fj, and 4-t- 
Thus, there is the entropy fee In 2 associated with the 
IPB. 

In general case we may introduce the free energy of the 
IPB as 



Gh — Hh — TSh 



(21) 



where is the enthalpy of the boundary creation and 
Sh its entropy. The IPBs break long range correlations 
between different parts of the system. Therefore, the 
correlations extend at the distances which are inversely 
proportional to the IPB concentration. The latter can be 
estimated &s~2- 



c b = e-^ 



(22) 



As can be seen from the explicit ID solution above, c b 
also defines the width of the region around p c where the 
fast change of the adsorption isotherm takes place. This 
can be visualized with the use of the variable x introduced 



as 



/3(p - p c ) = c b x. 



(23) 



With the use of this variable one can establish on the ba- 
sus of Eq. (|17p the explicit form of the universal function 
in Eqs. @ and © as 



W(x) 



(24) 



We note, that W(±oo) = ±1, as necessary. From practi- 
cal point of view, Eqs. (0) and (fT?) are not very conve- 
nient for universality checks because both sides of these 
equations turn into zero at the critical point p = p c . This 
means that in the data measured or calculated with fi- 
nite precision on a non-singular background the universal 
behavior can be obscured by the errors. The singular be- 
havior can be considerably enhanced by differentiation 
with respect to the gas pressure (see Fig. 1 in Ref. H3) or 
with respect to the variable x: 



dp 
dx 



1 



2(1 + x 2 ) 3 / 2 ' 



(25) 



As will be shown in Sec. IIV1 this expression indeed de- 
scribes the isothermal compressibility of coherent de- 
posits on SWNTs near the steps of the adsorption 
isotherms. 



III. THE MODEL 

The configuration of a coherent deposit consisting of 
identical atoms or molecules in the submonolayer cover- 
age regime can be fully characterized by the occupation 
numbers ni — 0, 1 of the deposition sites i = 1,N. The 
configuration energy of the deposit can be expanded into 
an infinite series of effective cluster interactions (ECIs) 



,15-17 



= (E ads -p)J2 n * + J2 v a 



l>3 



+ V ijk nin i nk 

i>j>k 



(26) 



where we assumed that the system is homogeneous so the 
adsorption energy E ac i s is the same at each site; also, as in 
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FIG. 1: Effective cluster interactions used to fit Hamiltonian 
in Eq. (|26[) to the ab initio data of Ref. @- The triangular lat- 
tice of deposition sites is shown (the deposition site is defined 
as the center of the carbon hexagon, see Fig. [2}. The tube 
axis is directed vertically. The leftmost and the rightmost 
sites in the odd rows on the drawing represent the same site 
on the tube. 



Eq. (|TU)) . we included into the Hamiltonian the chemical 
potential /i to control the coverage. To find the inter- 
action parameters in Eq. (|26[) . one needs, according to 
established methodology^^— to compute the energies of 
a sufficiently large number of different adsorbate struc- 
tures and then fit these energies to the lattice gas Hamil- 
tonian (|26p with sufficient number of ECIs. The energies 
are usually calculated ab initio but model calculations 
present viable alternative ) 15-17 Yet another possibility is 
to adjust the interactions to the experimental data (see 
the discussion and the bibliography in Ref. [Tot) £2r— 

In our calculations below we consider the adsorption of 
potassium on the surface of the zig-zag (6,0) carbon nan- 
otube studied in Ref. [H where the energies of six ordered 
structures were calculated within an ab initio approach. 
We remind that potassium and other metal deposits are 
directly related to the problem of hydrogen storage4~— ^ 

From the structures considered in Ref. |fj one can de- 
duce that at least third neighbor pair interactions need be 
included into the Hamiltonian (|2"6"|) (see their Fig. 1 and 
our Figs. [T]-[5]). Because of the tube anisotropy, the num- 
ber of pair interactions is, in fact equal to six, as shown 
in Fig. [TJ This is equal to the number of energy values 
we have at our disposition which is insufficient even to fit 
the pair interactions because we need also to determine 
the adsorption energy E a d s - 

To overcome this difficulty we, following Ref. |H, as- 
sume that the pair interactions between the potassium 
atoms can be approximately described by the Morse po- 
tential 



V(r l3 ) = e ( e -2«(n 3 -r ) _ 2e- a ( r *i- r ti) 



(27) 



which depends on three parameters e, a, and r$. Tak- 
ing into account E a d s , we are left with the possibility to 
adjust two more parameters. This turned out to be indis- 
pensable because the data of Ref. could not be fitted to 
the Hamiltonian containing only pair interactions. This 
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FIG. 2: The ground state structures on the surface of (6,0) 
SWNT with fillings 1/4, 5/12, and 2/3 found in Monte Carlo 
simulated annealing described in the text. The small dots 
correspond to the carbon atoms and the large dots the potas- 
sium atoms. The ranges of stability of these structures with 
the change of the chemical potential are shown on Fig. [3] 
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FIG. 3: Zero temperature phase diagram of potassium ad- 
sorbed on the surface of (6,0) carbon nanotube derived on 
the basis of Hamiltonian in Eq. (|26[) fitted to the data of Ref. 
@. The dotted lines represent the fit to some of the struc- 
tures found in that reference. The solid line corresponds to 
the new phase shown in Fig. [2] Ef is the energy of formation 
of the adsorbed structure and uk the chemical potential of 
potassium. The vertical line corresponds to bulk potassium. 6 



can be shown by expressing the energies of the structures 
in terms of (unknown) pair interactions and then estab- 
lishing exact relations between some of the energies by 
excluding the pair interactions from the expressions. The 
sum rules obtained in this way are strongly violated by 
the data of Ref. [f| 

Therefore, following Refs. we added two trio 

interactions comprising closely spaced atomsj 17 ' 46 i 47 By 
trial and error procedure we were able to achieve a very 
accurate fit to the six energies with two trio interactions 
shown in Fig. [T] and with the parameters presented in 
Table H 

The pair interactions presented in the table were ob- 
tained with the following parameters of the Morse po- 
tential: e = 0.136 eV, r = 5.69 A, and a = 0.426. This 
can be compared with the values obtained for the adsorp- 



TABLE I: Interactions entering Hamiltonian in Eq. (|26[) (eV) 



(/ 


V? 


vi 


vi 




Eads 


a 


-0.1072 


-4.119-10" 2 


-7.49-10" 2 


0.242 


-1.256 


b 


0.3427 


-0.1084 


-0.1268 


-2.825-10" 2 
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tion of potassium on copper: 43 e = 0.466 eV, r = 6 A, 
and a = 0.66. Taking into account that the systems are 
very different, our estimates look reasonable. A some- 
what too small value of e which define the attractive in- 
teraction between the potassium atoms may be due to 
the Coulomb repulsion because of the considerable and 
strongly coverage dependent charge transfer between the 
potassium and the substrate i 19 i 21 The value of the ad- 
sorption energy in Table U also agrees well with recent ab 
initio estimates^ 

Because our statistical approach is based on the grand 
ensemble, ECIs in Eq. ([26| do not depend on the coverage 
p which does not enter as a parameter in the formalism 
but is a dependent quantity calculated according to Eq. 
([5]). The concentration-independence may look unphysi- 
cal because the charge transfer which strongly influences 
the Coulomb interatomic interaction strongly depends on 
coverage.— Besides, in a similar problem in binary alloys 
it was shown that in the canonical formalism ECIs do de- 
pend on the concentration* 1 ^^ In Refs. |48||49| . however, 
it was shown that, if properly implemented, both for- 
malisms are equivalent. Formally in the grand ensemble 
the concentration independent cluster interactions (the 
pair ones, the three-body and higher) cooperate to re- 
produce the concentration dependence of the pair inter- 
actions of the canonical formalismi 48 i 49 

Physically the need for the three-body and higher ECIs 
can be understood as follows. In the case of only pair 
interactions there exists a "particle-hole" symmetry 



-ffpair = Vi^njnj-n rij = y~] Vijhjhj-fi' ^ hj+C, 

ij i i>j i 

(28) 

where hi = 1 — rij, // is a renormalized chemical poten- 
tial, and C a configuration- independent constant. Be- 
cause the chemical potential is an adjustable parameter 
fixing the coverage, in the case of constant pair interac- 
tions Vij it follows from Eq. (|28l) that the free energies 
calculated for coverage p and 1 — p differ only by the 
constant C. Thus, the derivatives of the free energy with 
respect to p whose singularities correspond to phase tran- 
sitions (at least, in 2- and 3D systems) are distributed 
symmetrically with respect to p = 1/2. This means that 
the phase diagram of the system with only pair interac- 
tions is strictly symmetric* 1 ^ But physically this is rarely 
the case, so the presence of higher ECIs is very common. 
Quite often the asymmetry of the diagram is strong which 
require the presence of large three-body ECIs compara- 
ble in magnitude with the pair interactions. 14 As can be 
concluded from the value of the trio interaction V^ 1 in 
Table HI this is also the case in the potassium adsorbates 
under consideration. 



IV. POTASSIUM ADSORPTION ON THE (6,0) 
SWNT 

The model of potassium adsorption on the (6,0) SWNT 
considered in previous section can be solved with the use 
of the same TM technique as in Sec. Ill Al only the TM 
will be much more complex than Eq. (1131) . To account 
for all interactions shown in Fig. [T] we need the TM of a 
rather large size 2 14 = 16384, as explained in Appendix. 
Fortunately, only the largest eigenvalue is needed for our 
purposes so the efficient technique of finding extremal 
eigenvalues of nonsymmetric matrices due to Arnoldi as 
realized in the software package ARPACK 50 could be used. 
Defining the reduced (per site) free energy 

<f> = -/3- 1 \n'Z/N, (29) 

the coverage can be calculated as [see Eq. (jTTJ)] 

p = J2( n i)/ N = - d( l , /d»- (30) 

i 

The adsorption isotherms for three different tempera- 
tures shown in Fig. 2] were calculated according to this 
definition with the use of the Hellmann-Feynman theo- 
rem to improve precision (see Appendix). As noted ear- 
lier, in the calculations we used Hamiltonian (|26[) with 
parameters from Table HI Though the parameters fitted 
the data of Ref. very accurately, in our calculations we 
did not see the quasi-transitions at or close to the val- 
ues of the chemical potential shown at Fig. 2 of Ref. 0. 
Our TM solution, however, is exact up to the computa- 
tional errors. Therefore, to establish the source of the 
discrepancy, we took the values of coverages at the low- 
est temperature (400 K) curve in Fig. [5] which to a high 
accuracy were equal to 1/4, 5/12, and 2/3 and performed 
Monte Carlo simulations in the framework of the canon- 
ical ensemble at these coverages.— The simulations with 
the use of the Metropolis algorithm were started at high 
temperatures and the system was gradually annealed to 
its ground state. At the coverages 1/4 and 2/3 we recov- 
ered the structures of Ref. while at coverage 5/12 an 
additional structure shown in Fig. [5] was found. It turned 
out to have lower energy than their structure at p = 1/2. 
This is shown in our Fig. |3]which is to be compared with 
Fig. 2(a) of Ref. S 

From the point of view of the CEM, the appearance of 
a ground state unaccounted for in ab initio calculations 
diminishes the accuracy of the whole scheme because the 
ground states are the only ones directly observable in 
statistical calculations (as temperature tends to zero)^ 
Therefore, it is highly desirable that they entered into 
the set of the structures calculated ab initio. While this 
point is important for the accuracy of the approach, in 
the present paper our main interest is in the universal 
features of the thermodynamics which do not depend on 
the accuracy of the Hamiltonian. So we believe that as 
long as the order of magnitude of the interactions are as- 
sessed correctly, the parameters of Table Q] are sufficiently 
adequate for our purposes. 
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FIG. 4: Adsorption isotherms at different temperatures for 
the system described by the Hamiltonian in Eq. (|26[1 with 
parameters given in Table [I] Upper curves are shifted up by 
0.3 with respect to the preceding curve for better visibility. 



Thus, according to our fit the (2 x 2)i?0° structure from 
Fig. [21 is the ground state of Hamiltonian (|2l>| in the in- 
terval of the chemical potentials -1.585 eV < px < -1-541 
eV (see Fig- EJ) j the 5/12 structure is stable for -1.541 eV 
< n K < -1.243 eV and the (a/3 x \/3)i?30° structure with 
p = 2/3 is the enrgy minimum for p^ larger than -1.243 
eV and up to the chemical potential of the bulk potas- 
sium calculated in Ref. H to be equal to -1.15 eV (the 
vertical line in our Fig. [3]). At finite temperature the 
zero-temperature boundaries between the ordered struc- 
tures give rise to three quasi-transition steps seen in Fig. 
UJ For simplicity we will refer to these transitions in order 
of their appearance from left to right as (quasi)transition 
number one, two, and three, respectively. 



A. Isothermal susceptibility and the universality 

The expression for the susceptibility with respect to 
the change of the chemical potential 



dp/dp = P^2((rii - p)(n - p)) 



(31) 



can be derived form Eqs. (| 10[) and (fTTj) . It can be used 
to assess the correlation length which we need in Eq. . 
As is known, both the susceptibility and the correlation 
length diverge at critical pointSi 26 i 31 ' 32 Thus, the points 
of the quasi-phase transitions in ID at finite temperature 
can be identified as the maxima of the correlation length, 
as suggested in Ref. El Below we will determine in this 
way the value p c of the critical chemical potential. Be- 
sides, Eq. (|3"Tj) can also be directly related to the isother- 
mal compressibility because at constant temperature the 
chemical potential is proportional to the logarithm of the 
pressure of the ideal gas. 12 Furthermore, with the use of 
Eq. ([7]) this quantity can be directly connected with the 
magnetic susceptibility of the Ising model. 

As can be seen from Fig. [U the quasi-transitions at 
the lowest temperature are so steep that can be easily 



0.01 r 




0.001 r 



FIG. 5: Low temperature behavior of the susceptibility Eq. 
(|31[) during three quasi-transitions seen in Fig. [4] The param- 
eters defining are presented in Table [XT] Ap is the change 
of the coverage during the transition which is 1/4 in the first 
and the third transition and 1/6 in the second one. The data 
plotted were calculated at temperatures T = 406 K (♦), T = 
464 K (A), and T = 580 K (o) near the first transition point, 
at T = 406 K (x) and T — 580 K (□) near the second tran- 
sition point, and at T = 464 K (0) and T = 580 K (+) near 
the third transition point. 



confounded with the true first order transitions. A pos- 
sible check on whether the transition is the true one is 
via the susceptibility in Eq. (j5T[) which should diverge at 
the true phase transition. We calculated this quantity by 
numerical differentiation of p{p). The results are plotted 
in the form of Eq. (|2"5"|) in Fig. [SJ The parameters corre- 
sponding to the IPBs are presented in Table HH As can 
be seen, all values in the table are reasonable from the 
point of view of the IPB interpretation: the entropies are 
all greater than the lower bound Ub In 2 rj 0.69/cb of the 
purely ID model of Sec. Ill Al and all enthalpies are no- 
tably larger than the individual interatomic interaction 
energies in Table HI We, however, were unable to calcu- 
late these values on the basis of an IPB model. Qual- 
itatively it is clear that IPBs in our system correspond 
to the rearrangement of the atoms from one phase to 
another and taking into account the complexity of some 
of them (see our Fig. [2] and Ref. @) the boundary may 
be not easy to guess. But even in the simplest case of 
the quasi-transition 1 between the empty lattice and the 
(2 x 2)i?0° phase characterized by six third-neighbor cou- 
plings between the successive atomic layers (see Fig. [5]), 
the naive calculation Eb = 61^/2 [by analogy with Eq. 
(|2TJ)) ] gives only 0.38 eV instead of 0.435. Also the entropy 
2.4&B suggests that the interface is rough, as is usual in 
2D Ising-like systems^ In our Monte Carlo simulations, 
however, the IPBs were very flat, at least at low temper- 
atures. Below we discuss some other possibilities which 
would require, however, additional investigations of this 
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TABLE II: Fitted values of the enthalpy and the entropy en- 
tering Gb in Eq. (|2ip for the three quasi-transitions seen in 
Fig.g] 



Quasi-transition No. 


H b (eV) 


Sb{k B ) 


1 


0.435 


2.40 


2 


0.545 


4.34 


3 


0.420 


0.80 



V. DISCUSSION 

In the present paper with the use of numerically ac- 
curate technique we were able to resolve the very steep 
behavior seen in the isotherms of adsorption on the nan- 
otube surfaces at low temperatures. 8 Such behavior is of 
considerable interest both from practical and from funda- 
mental points of view. On the one hand, it describes the 
response of the system to small variations in the exter- 
nal parameters; on the other hand, the strong response 
to the changes of the parameters can provide accurate 
information about microscopic interactions. 

In our calculations we used a realistic lattice gas model 
containing six anisotropic pair interactions and two clus- 
ter interactions among atomic trios derived in the frame- 
work of the the cluster expansion technique 1 ^— on the 
basis of ab initio electronic structure calculations.^ 

Despite the complexity of the model, its critical be- 
havior turned out to be the same as in the ID Ising (or, 
equivalently, lattice gas) model with nearest neighbor in- 
teractions. This agrees with the universality hypothesis 
of the theory of critical phenomena yet is a non-trivial re- 
sult because contrary to 3D case ) 29 i 30 in ID this hypoth- 
esis cannot be justified in the framework of the renormal- 
ization group approach for Hamiltonians with arbitrary 
interactions.— A qualitative explanation may be based 
on the very long range correlations present in the sys- 
tem at low temperature. The correlation length which 
can be assessed from the right hand side of Eq. (|3 1 1) is 
of 0(1/ Cb) and reaches values of O(10 4 ) as can be esti- 
mated from Table [III This means that the structures are 
correlated at very long distances and, using the language 
of the renormalization group and the Ising model, the 
block spin transformation can be efficient in theoretical 
description of the system. Because the tube diameter is 
much smaller than the correlation length, the block spins 
will comprise all the spins around the tube circumference 
as well as considerable block of sites along the tube. In 
this picture the interactions of sufficiently short range 
will connect only the nearest neighbor block spins thus 
making the system effectively equivalent to strictly ID 
model with only nearest neighbor interactions. 

The correlation length in O(10 4 ) of lattice spacings 
along the tube means that the whole nanotube can be 
covered by the ordered structure at temperatures as high 
as 400 K. (We note that only in infinite systems the long 
range order should be broken in lDj 1 ^ in a finite sys- 
tem the order can extend along the whole tube length^) 



Thus, from the hydrogen storage standpoint, at ambient 
conditions the potassium structure can be treated as an 
inert (in statistical sense) substrate while the hydrogen 
molecules treated within a statistical approach which can 
be based on the formalism developed in the present pa- 
per. 

In this study we concentrated on the universality for 
the following reasons. First, because we had at our dispo- 
sition the energies of only six ab initio calculated struc- 
tures, the accuracy of the cluster Hamiltonian was rather 
poor. Therefore, only orders of magnitude of the quan- 
tities of interest could be calculated judging from the 
fact that even 60 structures calculated in Ref. [ID did 
not allow to calculate the phase transition temperatures 
with accuracy better than 50%. The universal behavior, 
however, is the same for all Hamiltonians belonging to 
the same universality class. The second reason was that 
in many cases the surface structures are not commensu- 
rate with the substrate Yet they can be as good hy- 
drogen adsorbers as the commensurate structures. But, 
as we noted in the Introduction, the lattice structure 
is not needed for the critical behavior to be universal, 
as the gas-liquid transitions in 1-, 2-, and 3D systems 
showier— According to Ref. HH the liquid can be viewed 
as a crystalline state filled with topological defects, such 
as dislocations and disclinations. The same can be said 
about the incommensurate surface layers^ Therefore, 
one might expect that the universal critical behavior may 
take place also in the incommensurate cases. This predic- 
tion should be amenable to experimental verification on 
the isotherms of SWNTs covered with incommensurate 
phases of inert gases4 

The third reason for studying the universality was that 
while the TM technique is an accurate and efficient tool 
for treating the ordering of coherent adsorbates on sur- 
faces of small nanotubes, the size of TM grows exponen- 
tially with the tube diameter and with the range of inter- 
actions. This means that for only slightly larger tube or 
longer-ranged interaction the TM will became unmanage- 
ably large. There exist viable alternatives to the TM in 
solution of this kind of problems: the mean field approxi- 
mation and especially the Monte Carlo methodi 14 ' 52 ' 57 i 58 
Both techniques, however, meet with difficulties in treat- 
ing the fine details of the abrupt phase transition-like 
changes seen on the adsorption isotherms. The results 
obtained in the present paper are aimed at resolving this 
difficulty. The mean field or the Monte Carlo methods 
can accurately predict the position of the transition while 
the universality in the transition curves observed in our 
study should provide its fine details. There remains the 
problem of finding the values of the two parameters which 
describe the behavior quantitatively. A block-spin renor- 
malization group and/or the low-temperature expansion 
are probable candidate tools for attacking this problem. 
Further work is needed to clarify this point. 
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function. This is because when adding a site to the sys- 
tem consisting of K > R sites only the interactions with 
the last R sites need be taken into account. The ac- 
counting can be done with the use of the vector partition 
function Z^ K > whose components are the partial traces 
over all except the last R sites (the sites are numbered 
from right to left) 



Appendix: Sparse transfer matrices 

Our TM approach belongs to the general category of 
TM methods based on sparse matrices initiated in Refs. 
I35ll36l : further bibliography can be found in Ref. H^. The 
problem of adsorption on triangular lattice with account 
of the second neighbor and trio interactions was previ- 
ously studied within similar framework in Ref. [3 but no 
details were given. We believe that our technique pre- 
sented below is particularly simple and easy to use. 

The advantage of using sparse TMs is that instead of 
I 2 matrix elements of a conventional dense I x I TM ma- 
trix (see, e. g., Ref. [T3) one deals with matrices contain- 
ing only 0(1) nontrivial entries. Because the size of TM 
scales with the range of interactions R exponentially as 



(A.l) 



and in practical calculations reaches significant values (e. 
g., 2 14 = 16384 in the present study), the gain in numer- 
ical efficiency from using sparse TMs can be enormous. 

The interaction range R in Eq. (|A.1[) for Hamiltonian 
([25)1 is defined as the longest range of the cluster inter- 
actions it contains. The range of a cluster interaction 
Vj. ,■ rii. ... rii is defined as 



R = L 



(A.2) 



where i ma x and i m in are the maximal and the minimal 
indices amonf i±, . . . ,i m . For example, the cluster inter- 
action nini+ifii+2 has the range R = 2. 

The finite range of interactions in the Hamiltonian 
makes possible a recursive calculation of the partition 



Z( K ) = Tr 



which can be visualized as 



(A.3) 



(A.4) 



where the empty and filled circles correspond to the 
empty (rii = 0) or filled (rii = 1) sites in Eq. (|A.3[) 
while asterisks denote the sites over which the trace over 
the two possible values of filling has been taken; H^ K > 
is Hamiltonian (f2l)|) for a if-site system. The partition 
function is found from (|A.3I) as 




(A.5) 



Here the bar over the number denotes that its binary 
representation is meant 



A = (a R ^ 



. aia )j 



(A, 



where = 0, 1 correspond to the filling of site k. The 
subscript R reminds that the term within parentheses is 
the binary representation, not the product, and that its 
length is equal to R. For example, T = 00 ... 01 with 
R — 1 zeros before the unity means that there is R — 1 
empty sites before the filled one. 

The general form of the TM can be understood from 
the recurrence equation 
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(A.7) 



where the column vectors correspond to Z^ N 1 - 1 and with the components denoted by their subscripts in Eq. 
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(|A.4I) for brevity. The subscript N of the TM is the site 
index for all b a entering the matrix. We note that we 
use the same symbol N for the system size and for the 
recurrent relation in order to stress that at every itera- 
tion we obtain the (vector) partition function of a system 
of size N, i. e., that the partition functions obtained at 
intermediary steps are not in any way deficient. 

The structure of TM in (IA.7I) is physically transparent. 
Having added site N to the system consisting of N — 1 
sites we first have to account for the interaction of this 
site with the rest of the system and then take the trace 
over the (N — R)-th site because with the radius of inter- 
actions being R all interactions of this site with the rest 
of the system have already been taken into account. Tak- 
ing the trace amounts to adding with appropriate weights 
two Z^ N ~ X > differing by the filling of site N — R. In the 
case when site N is empty the weights are equal to unity 
because the empty site does not interact with anything 
and the interaction energy is zero. These terms occupy 
the upper half of the TM (|A~7)) . The lower half of the 
matrix contains the terms corresponding to the interac- 
tion of the occupied site N with the rest of the system. 
The term 

b aN = exp(-/3A£ aA r) (A.8) 

is the Boltzmann weight corresponding to the interaction 
of the atom at site N with the configuration of atoms 
corresponding to Z^~ 1 '] AE & n in (|A.8[) is the energy 
of interaction of the atom at site N with configuration a 
on sites N - 1, N - 2, . . . , N - R. 



1. Application to adsorption on (6,0) nanotube 

In Fig. [5] are shown both the enumeration of sites we 
chose for the (6,0) nanotube and the fourteen sites along 
the path with which atom at site i+17 can interact if they 
are also filled with atoms. The furthest neighbor site i + 3 
is defined by the the 3rd neighbor interaction V£ which 
can reach it (see Fig. [TJ. Thus, according to Eq. (|A.2[) 
R = 14 and the size of our TMs is 2 14 = 16384. This is 
the number of configurations we need to account for in 
our transfer matrices. From Fig. [5] one can see that as 
the sites are being added one after another in the top row 
the relative placement of the 14th neighbor change with 
respect to the added site. Because of this the transfer 
matrices for neighbor sites are different, except for sites 
i + 15 and i + 16 which is due to the particular interac- 
tions entering our Hamiltonian. This can be seen from 
Table ITO1 where the pair interactions accounted for in the 
Boltzmann factors bk entering the TMs are presented. 
Similar tables can be composed for the trio interactions. 

It is easy to see that the structure of the TMs repeats 
after each six steps, for example, when the row gets filled. 
This allows one to compute the reduced free energy of 
the system Eq. (|T4"|) through the logarithm of the largest 



i + 13 « + 14 « + 15 i + 16 i + 17 « + 18 (i + 13) 




i+1 i + 2 i + 3 i + i i + o i + 6 (i + 1) 



FIG. 6: Enumeration of sites on the (6,0) nanotube used in 
the construction of the TM. Black dots denote the fourteen 
neighbors of site i + 17. The 14-th neighbor i + 3 is defined 
by the interaction V3 in Fig. \T\ 

TABLE III: Interactions of atoms in the top row on Fig.[6]with 
atoms on fourteen preceding sites (index i has been omitted 
for brevity). The arrows point to the value they represent. 



Neighbor No. 
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14 





«- 


vi 
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<— 



eigenvalue A+ of the product of six TMs, e. g., of those 
corresponding to sites from i + 13toi + 18in Fig. [5] as 

(p = -fc B TlnA + /6. (A.9) 

2. Adsorption isotherms 

To draw the adsorption isotherm one has to calculate 
the derivative of <j) with respect to the chemical poten- 
tial [see Eq. (|3"0"|) ]. With the fast variation of the deriva- 
tive in the most interesting region in the vicinity of the 
quasi-transition, the numerical differentiation can be un- 
reliable. More accurate results can be obtained with the 
use of the Hellmann-Feynman theorem: 

# = 1 (+\dMe/d(M\+) (A1Q) 
P d^i 6 (+|M 6 |+) 

where Mq is the product of the six TMs, as explained 
above, |+) is the eigenvector corresponding to A+, and 
(+| the eigenvector of the transposed matrix because our 
TMs are not symmetric. Due to the simplicity of our TMs 
the derivative in Eq. (|A.10[) is very easy to calculate: the 
upper part of each of the six TMs entering Mq should 
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simply be successively set to zero while the lower part change the exponential function exp(/3/x). 
remains the same because the differentiation docs not 
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